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Abstract 

This  paper  documents  the  construction  of  a finite  element,  called  the  Sibson  element.  The  shape 
function  of  this  element  is  formed  on  rectangular  grids  by  C1  splines  defined  on  a triangulation  of  each 
subrectangle  by  dividing  it  into  four  subtriangles  formed  by  drawing  the  diagonals.  The  splines  are 
constructed  from  bivariate  cubic  polynomials  z(x,  y)  and  are  written  in  such  a way  that  they  are  linear 
functions  of  the  values  2,  dzjdx , dz/dy  at  each  node  of  the  rectangle  with  bivariate  polynomial  coeffi- 
cients up  to  order  three.  Conditions  are  given  for  the  existence  of  such  an  element.  They  are  used  to 
construct  the  bivariate  polynomial  coefficients,  first  for  a unit  rectangle  and  then  for  a general  rectangle. 
Since  the  first  and  second  derivatives  of  these  functions  are  sometimes  needed  they  are  also  given. 

Key  Words:  bivariate  polynomial;  rectangular  mesh;  shape  function;  Sibson  element;  spline;  triangu- 
lation 


1 Introduction 

When  modeling  urban  and  rural  terrain  data  we  are  usually  given  points  {xz,yl,zz)  for  i — 1,2,-- -,  not 
necessarily  uniformly  spaced.  We  would  like  to  construct  a C1  function  to  interpolate  the  given  data. 
Lawson  [5,  6]  has  pioneered  C1  interpolation  over  triangular  grids.  Han  and  Schumaker  [1]  interpolate 
data  on  rectangular  grids  by  C 1 splines  defined  on  a triangulation  of  each  subrectangle  by  dividing  it  into 
four  subtriangles  formed  by  drawing  the  diagonals.  This  subdivision  is  called  the  Sibson  split.  They  then 
construct  the  C 1 spline  by  using  a Bernstein-Bezier  form.  Recently  Lavery  [3,  4]  has  used  the  Sibson  split 
on  rectangular  grids  of  terrain  data  but  has  constructed  Cl  splines  over  these  Sibson  elements  directly  from 
bivariate  cubic  polynomials.  He  briefly  describes  the  derivation  in  Lavery  [4]  but  only  gives  the  final  formula 
for  the  spline  on  one  of  the  subtriangles  of  a subrectangle.  It  is  the  object  of  this  report  to  present  the  details 
involved  with  the  derivation  of  these  bivariate  cubic  elements. 

We  begin  with  a tensor-product  grid  with  nodes  xz,  for  i = 0, 1,  • • • , 7,  and  yJZ  for  j — 0, 1,  • • • , J,  that 
form  strictly  monotonic  partitions  of  the  finite  real  intervals  [xo,x/]  and  [.Vo,yj]  respectively.  The  nodes 
need  not  be  uniformly  spaced.  The  cubic  splines  z(x,y)  will  consist  of  C1  smooth,  piecewise  cubic  Sibson 
elements  on  the  tensor-product  grid.  We  also  suppose  that  we  are  given,  by  some  means,  elevation  and 
derivative  values  Zy,  z * = dzjdx,  zjf  = dz/dy  at  each  node. 

To  create  a Sibson  element  on  a rectangle  (xj,x,+j)  x (t/7,  yJ+1)  we  first  divide  the  rectangle  into  four 
triangles  by  drawing  the  two  diagonals  of  the  rectangle.  The  Sibson  element  is  a shape  function,  z(x,y), 
that  is  cubic  within  each  triangle,  is  C1  on  the  lines  separating  the  four  triangles,  has  dz/dx  that  is  linear 
along  the  edges  x = xz  and  x — xz+\  of  the  rectangle  and  has  derivative  dz/dy  that  is  linear  along  the  edges 
V = V]  and  y = yJ  + 1. 
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The  report  is  structured  as  follows.  Section  2 describes  the  specific  sufficient  conditions  for  constructing 
a Sibson  element  on  a square  with  unit  length  on  each  side.  The  necessary  equations  for  determining  the 
bivariate  polynomials  in  each  of  the  four  Sibson  triangles  of  the  unit  element  are  derived  in  Section  3.  Solving 
for  the  bivariate  polynomial  equations  is  done  in  Section  4.  In  Section  5 the  coeffficients  are  written  in  terms 
of  the  three  items  of  information  at  each  of  the  corner  nodes.  In  Section  6 the  final  Sibson  element  shape 
functions  are  constructed  for  the  unit  element.  In  Section  7 the  Sibson  element  shape  functions  on  the  unit 
square  are  transformed  to  shape  functions  on  an  arbitrary  scaled  rectangle.  Finally,  in  Section  8 the  first 
and  second  derivatives  of  the  Sibson  element  splines  are  computed. 

2 Defining  the  Conditions  for  Constructing  a Unit  Sibson  Element 

In  this  section  we  develop  the  Sibson  element  on  a reference  unit  square  in  order  to  simplify  the  construction. 
In  a later  section  we  will  introduce  a transformation  that  will  map  this  element  to  non-uniform  Sibson 
elements.  Figure  1 shows  a unit  square  with  positive  r-axis  pointing  downwards  and  positive  s-axis  pointing 
to  the  right.  The  corners  are  identified  as  (0,0),  (1,0),  (1, 1),  (0,  1).  The  diagonal  line  from  (0,0)  to  (1, 1)  is 
given  by  the  equation  s = r and  that  from  (1,0)  to  (0, 1)  is  given  by  s = r — 1.  The  triangles  are  labeled 
0, 1,2,3,  beginning  with  the  triangle  identified  by  the  points  (0, 0),  (1,  0),  (1/2, 1/2),  and  proceeding  in  a 
counterclockwise  manner. 


s 


Figure  1:  Unit  Sibson  Element 

Cubic  polynomials  are  defined  on  each  triangle  in  the  following  two  forms.  First,  for  triangles  indexed 
by  k = 0,  ■ • • , 3 the  bivariate  polynomials  are  written  as 

zk(r,  s ) = c£0r3  + c%xr2s  + cx2rs 2 + c£3s3 

+ c207'2  + cllr5  + c02s2  + c10r  (!) 


2 


+c01s  + c00 


and  for  triangles  indexed  by  k = 1,  2 they  are  written  as 


Zk(r,s)  - 


(1  - r) 


l30 

+c2o(l  ~ rY  + cn l 


+ c2i(l—r)  (1  — s)  + C12(1  — r)(l  ~~  s)  +c03( 


1-5) 


1 - r)(l  - s)  4-  cg2(  1 - s)2  + cf0(l  - r) 


(2) 


+coi  (1  — s)  + c. 


oo- 


These  bivariate  polynomial  forms  are  selected  so  that  a point  (r,  s)  in  triangle  0 is  matched  with  a 
symmetric  point  ((1  — r),  (1  — s))  in  triangle  2 along  a line  through  the  point  (1/2,  1/2).  A similar  symmetry 
holds  for  triangles  1 and  3.  The  coefficients  are  indexed  as  follows.  The  superscript  k indexes  the  triangle 


and  the  subscript  pq  represents  the  exponentials  in  the  term  Cpqrpsq  or  c£q(  1 — r)p(  1 — s)q . 

These  four  bivariate  cubic  polynomials  must  interpolate  at  the  four  corners,  must  have  ^ that  is  linear 
along  the  edges  x = and  x = xl+i  of  the  rectangle  and  have  derivative  dz/dy  that  is  linear  along  the 
edges  y = Vj  and  y = pj+\  and  satisfy  the  C1  conditions  across  the  diagonals.  The  interpolation  conditions 
at  the  corners  are  given  by  the  following  identities. 


2o(0,  0) 

dz0 


dr 


(0,0)  = 


23(0,0) 

dz3 


dr 


(0,0) 


tr(0’0) 

2o(l,0) 


2l(l,0) 

dzi 


dr 


(1,0) 


dz0 


(1.0)  = 


ds 

22(0,1) 
^(0,1)  = 
= 

2l(M)  = 

- 

- 


= 23(0,1) 

= >,D 

dr 


(3) 


dZ3 


(0,1) 


ds 

22(1,1) 

£(1.1> 


In  the  unit  cell  case,  the  outer  boundary  conditions  imply  that  the  following  partial  derivatives,  for  r € 

[0,1],  5 € [0,1], 

dz°  ( n\ 

a7(r'0) 


dzi 

M) 

dr 


(4) 


are  linear.  Finally,  the  polynomials  must  be  C1  across  the  diagonals.  Between  triangles  0 and  3,  along  the 
line  s — r,  the  following  equations  must  be  satisfied. 

z0{r,r)  = z3(r,r), 
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(5) 


dz0 

dr 

dz0 


, v dzz 

(r,r)  = ^r(r,r), 

dz3 


ds 

Between  triangles  1 and  2,  along  the  line  s = r, 


-^-(r,r)  = — (r,r). 


5s 


2j(r,r)  = z2(AH, 

dz2 


dz\  . 
•»(r’r)  = 


5zi 

5s 


(r,r)  = 


Between  triangles  0 and  1,  along  the  line  s = 1 — r, 


dr 

dz2 

ds 


r,  r), 
(r,r). 


zo(r,l-r)  = zj(r,l-r), 


dz0 


dz\ 


— (r,l-r)  = 

^(r,l-r)  = $i(r,l-r). 


ds 


ds 


Finally,  between  triangles  2 and  3,  along  the  line  s = 1 — r, 

z2(r,  1-r)  = z3(r,  1 - r), 

dz2  dz3 

— (r,l-r)  = -^(r,  1 - r), 

<922  , . v <9*3  , , v 

17(r'1_r)  = a7(r’1-r)' 


(6) 


(7) 


(8) 


All  of  these  conditions  will  produce  the  necessary  equations  that  must  be  satisfied  by  the  forty  unknown 
coefficients.  In  the  next  section  we  will  develop  these  equations. 


3 Necessary  Equations  for  Coefficient  Determination 

The  equations  developed  in  this  section  are  all  determined  by  the  conditions  in  the  previous  section.  They 
arise  from  simply  inserting  the  appropriate  polynomial,  (1)  or  (2),  into  either  the  interpolation,  boundary  or 
diagonal  conditions  and  equating  appropriate  coefficients.  In  this  section  zl0  represents  the  elevation  data  at 
the  corner  (z,  j)  of  a cell.  z[3  is  the  derivative  with  respect  to  r at  (i,  j)  and  z|  is  the  derivative  with  respect 
to  s at 

We  first  develop  the  equations  that  arise  from  the  interpolation  conditions.  We  will  use  (3)  to  develop 
the  coefficient  relations  that  arise  as  a result  of  the  interpolation  at  the  unit  cell  corners.  In  particular,  at 
(0, 0)  simple  substitutions  show  that 


r0 

c00 

= 2oo, 

(9) 

r3 

c00 

— 200, 

(10) 

r0 

c10 

II 

(11) 

c10 

II 

(12) 

C01 

— 200> 

(13) 

r3 

C01 

— 2oo- 

(14) 

At  (1,0)  the  coefficient  relations  are  given  by 

c30  + ^20  + 

c10  + c00  = 210 

(15) 
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At  (0, 1)  the  relations  are 


c03  + c02  + C01  + c00 

= 

z10 

(16) 

3c°  4-  9c°  4-  c° 
oc30  x-  ZC20  -h  C10 

= 

7r 

"10 

(17) 

_r  1 _ rl  _ rl 
c12  C11  c10 

- 

"10 

(18) 

ro  +r°  . ro 

°21  + C11  x-  C01 

= 

"10 

(19) 

— 3c1  — 2c1  — c1 
oc03  "^02  C01 

Z10 

(20) 

c30  + c20  + c10  + c00 

— 

Z01 

(21) 

c03  + c02  + C01  + c00 

= 

Z01 

(22) 

— 3clo  — 2c^0  — Cj0 

— 

Z01 

(23) 

4-  4-  r ^ 

c12  ^ C11  ^ c10 

= 

Z01 

(24) 

2 2 2 

°2\  ~~  C11  C01 

= 

s 

"01 

(25) 

3Cq3  + 2Cq2  + Cqj 

- 

"01 

(26) 

The  last  interpolation  relations,  given  at  (1, 1),  are 


ri 

coo  — 

Z11 

(27) 

r2 

c00  — 

211 

(28) 

rl 

c10  — 

~zu 

(29) 

r2 

cio  — 

~zn 

(30) 

C01  — 

~zii 

(31) 

r2 

C01  — 

~zu 

(32) 

At  the  external  boundaries  of  the  cells  dz/dx  is  linear  along  the  edges  x = xl  and  x = xl+1  of  the 
rectangle  and  the  derivative  dz/dy  is  linear  along  the  edges  y = yd  and  y = yJ+\.  By  differentiating  (1)  and 
(2),  evaluating  them  on  the  appropriate  boundaries  and  setting  the  quadratic  coefficients  to  zero  it  is  easy 
to  see  that 


ro 

= 0 

(33) 

rl 

c12 

= 0 

(34) 

4l 

= 0 

(35) 

c12 

- 0 

(36) 

The  C1  conditions  along  the  diagonals  are  obtained  by  simply  evaluating  the  relations  (5),  (6),  (7),  (8), 
and  equating  the  coefficients  of  like  powers.  Therefore,  along  s = r between  triangle  0 and  3 we  obtain  the 
equations 


r°  4_r°  4-  c°  4-c0 
c30  ^ ^ c12  ^ l03 

- 

C30  + C21  + c12  + C03 

(37) 

r°  J-,°  x,° 

*420  ' C11  x-  c02 

- 

4o  + C11  + c02 

(38) 

c°  + c° 
c10  ^ L01 

— 

c10  + C01 

(39) 

r0 

c00 

= 

r3 

c00 

(40) 

3^30  + 2c^j  + Cj2 

= 

3*^30  + 2c^!  + Cj2 

(41) 

2c^0  + Cji 

- 

2*^20  + C11 

(42) 

o 

O ^ 

o 

= 

c10 

(43) 

°21  + 2c12  + 3Cq3 

= 

C21  + 2Cj2  + 3Cq3 

(44) 

C1 1 + 2Co2 

= 

ch  + 2cq2 

(45) 

r0 

C01 

= 

r3 

C01 

(46) 
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Along  the  line  s — r,  between  triangles  1 and  2,  we  obtain 


C30  + 4l  + <02  + c03 
3^30  + 3(4j  + 3 c}2  + 3cq3 

— 

c30  + 4l  + Cj2  + Cq3 

(47) 

+ °20  + C11  + c02 

— 

3^0  + 3c21  + 3c22 
+ 3Cq3  + 4o  + C1 1 + 4)2 

(48) 

— 3c3q  — 3c2j  — 3c}2  — 3cq3  — 2c^0 

— 2c}  j — 2Cq2  — c}0  — Cqj 

= 

— 3C30  — 3c^j  — 3c  j2  — 3cq3 
— 2c|0  — 2Cjj  — 2Cq2  — Cj0  — Cqj 

(49) 

C30  + <41+  c12  + c03  + °20  + c\1 

+ c02  + c10  + C01  + C00 

= 

c30  + °21  + Cj2  + Cq3 
+ 4>0  + Cjj  + Cq2  + Cj0  + Cqj  + Cqq 

(50) 

_3c30  — 2c^j  — c}2 

= 

— 3c30  — 2c^j  — Cj2 

(51) 

6C30  + 4c21  + 2c}2  + 2c^0  + c}j 

= 

6C30  + 4c|j  + 2Cj2  + 2c|0  + Cjj 

(52) 

3c^q  — 2c21  — c}2  — 2 c^q  — cjj  — c}0 

= 

— 3c|0  — 2c2j  — Cj2  — 2c|o  — Cjj  — Cjo 

(53) 

— c21  — 2c}2  — 3cq3 

= 

~4l  — 2Cj2  — 3Cq3 

(54) 

2c^i  + 4c}2  + 6cq3  + c}j  + 2cq2 

= 

2c|l  + 4Cj2  + 6Cq3  + Cjj  + 2Cq2 

(55) 

~4\  — 2c}2  — 3Cq3  — c}j  — 2 Cq2  — 4l 

= 

~ C21  _ 2Cj2  — 3Cq3  — Cjj  — 2Cq2  — Cqj 

(56) 

Along  the  line  s — 1 — r,  between  triangles  0 and  1, 

we 

obtain 

C30  ~ C21  + Cj2  — C03 

-4)  + 4l  - cl2  + c03 

(57) 

4l  - 2Cj2  + 3Cq3  + C%0  — Cjj  + Cq2 

3c30  — 2<4j  + c]2  + 4 0 — c}j  + Cq2 

(58) 

c°  — 3r°  4-  r°  — 2c°  4-  c°  — c° 
c12  JC03  ^ C11  ~c02  ^ c10  C01 

- 

-3C30  + 4i  - 24o  + C11  - C10  + C01 

(59) 

r°  4 . r°  -A- r°  4- 
c03  c02  +c01+  c00 

= 

c30  + 4o  + 4o  + c00 

(60) 

3c30  — 2c21  + Cj  2 

= 

—3c30  + 2c^  j — c}2 

(61) 

2c§,  — 2cj2  + 2c2q  — Cjj 

= 

6C30  - 24\  + 24o  - c}j 

(62) 

c12  + C11  + c10 

= 

— 3c30  — 2 c^0  — c}0 

(63) 

°21  ~ 2Cj2  + 3Cq3 

= 

— 4i  + 2c}2  — 3cq3 

(64) 

2Cj2  — 6Cq3  + Cjj  — 2Cq2 

= 

24\  ~ 2ci2  4"  Cjj  — 2Cq2 

(65) 

3c03  4"  ^Cq2  + C01 

— 

4i  - cn  — coi 

(66) 

Finally,  along  the  line  s — 1 — r,  between  triangles  2 and  3,  we  obtain 

~c30  + °21  “ c12  + c03 

= 

C30  _ °21  4-  cf2  — Cq3 

(67) 

3c30  — 2c^j  + Cj2  + C^o  - Cjj  + Cq2 

= 

C21  — 2c12  + 3c03  + C20  — C11  + c02 

(68) 

“3c30  + C21  - 2 C^o  4-  Cjj  — Cj0  + Cqi 

= 

Cj2  — 3Cq3  + Cjj  — 2Cq2  + Cjo  — Cqj 

(69) 

c30  +<40  + c10  + c00 

— 

C03  + c02  4-  Cqj  + Cqo 

(70) 

— 3c30  + 2c^j  — Cj2 

= 

3c30  — 2c|j  + Cj2 

(71) 

6C30  _ 3c^j  + 2c^0  — Cjj 

= 

2c2j  — 2Cj2  — Cjj  + 2c^0 

(72) 

~3c|0  — 2c^0  — Cj0 

= 

Cl2  + C?J  + Cjo 

(73) 

— C21  + 2c12  — 3c03 

= 

C21  — 2Cj2  + 3c03 

(74) 

2cii  — 2cj2  + Cjj  — 2cq2 

= 

2Cj2  — 6C03  + Cjj  — 2Cq2 

(75) 

-c2  -c2  -c2 
*421  L11  C01 

= 

3Cq3  + 2Cq2  + Cqj 

(76) 
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4 Finding  a Solution 


There  are  sixty  eight  equations  developed  in  the  previous  section,  from  (9)  to  (76),  but  only  forty  unknowns, 
so  we  expect  twenty  eight  redundant  equations.  Of  the  forty  unknowns  sixteen  are  immediately  determined 
from  equations  (9)  through  (32)  and  (27)  through  (36).  With  these  given,  it  is  necessary  to  find  the  other 
twenty  four  solvable  equations  and  test  that  the  rest  of  the  equations  are  redundant.  Some  redundancies  are 
obvious.  For  example,  equation  (39)  is  redundant  with  (29)  through  (32),  whereas  other  redundancies  are 
not  so  clear. 

The  first  two  equations  that  can  be  solved  are  (15)  and  (17)  for  c20  and  c30  in  the  form 


ro  — 9„0  o„0 

c20  — a^lO  Z10  zc10  oc00 

r°  - r „ r°  r°  r°  r° 

c30  — 210  ~ 0*20  — c20  ~ c10  — c00 


(77) 


where  Cq0  and  Cj0  are  known  from  equations  (9)  and  (11)  respectively.  The  solutions  given  in  this  section 
will  be  developed  in  terms  of  coefficients  with  either  known  values  or  previously  computed  values.  The 
representation  of  the  coefficients  in  terms  of  the  data  at  the  four  corners  of  the  cell  will  be  done  in  the  next 
section.  That  step  of  the  process  required  the  assistance  of  a symbolic  manipulation  program. 

The  second  two  equations  that  can  be  solved  are  (16)  and  (20)  for  Cq2  and  Cq3  as 


c02  — + ^io  ~ 2coj  — 3cq0 

c03  = 210  — c02  — coi  — coo-  (78) 


The  next  two  equations  that  can  be  solved  are  (21)  and  (23)  for  c^q  and  c30  as 


-20 


— 32qi  T 2q,  2ciq  3c, 


'00 


10  Loo- 


c30  — 201  C20  c 

From  (22)  and  (20)  we  can  solve  for  Cq2  and  Cq3  as 

„3  or  ys  93  o„3 

c02  — a2oi  ^-oi  “c01  ^C, 


00 


„ „3 

3 

_3 

c03 

~01  — c02 

C01 

coo- 

ami  (24)  for 

r°  r 1 
cn,  cn, 

c\  j,  and  c3j 

n 

0 

0 

cn 

— z10 

C21 

C01 

A 

r 

A 

A 

cn 

— ~~z10 

~ c12 

~ c10 

2 

2 

2 

cn 

— — ^Ol 

°21 

~ C01 

„3 

r 

„3 

3 

C11 

— *01  " 

c12  ~ 

C10 

Since  c20,  Cjj,  and  c3,  have  been  computed  we  can  solve  (42)  for  c20  as 


c20  — c20  + o c 


2 11  2 

Again,  since  c3j,  Cq2,  and  c®j  have  been  computed  we  can  solve  (45)  for  Cq2  as 


1 

-c 


11  • 


r°  - Ir3  4-  r3 
c02  — 0C11  + c02 


2 11  2 

We  can  now  solve  (72)  for  c2 1 as 

G2i  = 9 (dc3o  — 2c21  + 2c20  — Cj,  + 2c32  + Cj,  — 2 c|0) 


(79) 


(80) 


(81) 


(82) 


(83) 


(84) 
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To  solve  for  c^q  insert  (51)  into  (53)  and  cancel  like  terms  in  (53)  to  get 

4o  = % (~cn  ~ c\o  + 24o  + cli  + cio)  (85) 

To  solve  for  Cq2,  multiply  (47)  by  3,  insert  the  result  into  (49)  and  cancel  like  terms  in  (49)  to  get 

c02  = 2 + 2cq2  + cj0  + Cqj  — 2c^0  — 2cjj  - c\Q  — c^)  (86) 

To  solve  for  Cj2,  first  solve  (66)  for  in  terms  of  known  values  and  the  unknown  Cq3,  which  will  be 
computed  later,  to  get 

c 21  = — C11  — C01  ~ 3c03  — 2Cq2  - Cq!-  (87) 

This  is  a temporary  equation  for  c^j  that  will  be  used  for  substitution  into  (64).  Once  this  substitution  is 
made,  terms  for  Cq3  cancel  in  equation  (64).  The  final  result  is 

ci2  — 2 (°2i  — cn  — coi  ~ 2cq2  — Cqj  — 2c]2  + 3 Cq3)  (88) 

where  all  of  the  terms  on  the  right  are  known. 

From  (63)  solve  for  c30  as 

c30  = g (~^c20  ~ cio  ~ ci2  ~ cii  ~ do)  (89) 

From  (61)  solve  for  c21  as 

4x  = \ (3c30  - 2c^!  + c\ 2 + 3c30  + c\2)  (90) 

From  (66)  we  can  now  solve  for  Cq3  as 

c03  = g (_c21  ~ cn  ~ co\  ~ 2cq2  — coi)  • (91) 

From  (41)  solve  for  c30  as 

c30  = 2 (3c30  + 2^21  + C12  — 2C2J  — Cj2)  . (92) 

From  (53)  solve  for  cf2  as 

c12  — 3c30  + 2c^j  + c\2  + 2<40  + cjj  + c|0 

— 3clo  — — 2c^0  — Cjj  — Cj0  (93) 

Finally,  from  (54)  we  can  solve  for  Cq3  as 

c03  = 2 + + 3c03  ~ ~ 2cj2)  • (94) 

As  is  clear  these  coefficients  were  solved  using  a subset  of  the  sixty  eight  equations.  All  of  the  other 
equations  are  redundant.  This  has  been  verified  by  a symbolic  manipulation  program. 
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5 Representing  the  Coefficients  in  Terms  of  Corner  Data 

The  coefficients  can  now  be  expanded  in  terms  of  the  twelve  data  items  at  the  corners  of  the  unit  square. 
This  again  is  done  by  a symbolic  manipulation  program.  The  results  are  given  below. 


c00 

■ 

z00 

c00 

= 

^00 

r0 

c10 

= 

~00 

r0 

C01 

= 

200 

c3 

c10 

= 

r 

Z00 

c3 

t01 

- 

s 

"00 

c00 

= 

2n 

A 

Loo 

- 

2ll 

A 

c10 

= 

~Z\\ 

A 

C01 

= 

~zh 

r2 

c10 

= 

~zii 

C01 

= 

-*'n 

4i 

= 

0 

A 

c12 

= 

0 

4 

= 

0 

c12 

= 

0 

A 

c20 

= 

3zio  — 2jq  — 2zqq  — 3^00 

ro 

c30 

= 

-2^io  + ^oo  + 220o  + zio 

A 

c02 

= 

3-io  + 2j0  + 2-Zij  — 3zn 

A 

c03 

= 

— 2sio  — -jo  — + 2—ii 

0 

c20 

= 

3-oi  + zoi  + 2zru  — 3zi  i 

r2 

c30 

-2z01-zr01-zru+2zu 

c02 

- 

3-oi  — zoi  ~~  2zqq  — 3zoo 

C03 

= 

— 2zoi  + 2qi  + Zoo  + 2zoo 

c°n 

= 

•Vs  _ -S 

"10  "00 

ch 

= 

T T 

Z11  — z10 

c2n 

= 

-s  _ 7s 

"11  "01 

c3n 

= 

r r 

"01  ~ "00 

c3 

c20 

= 

3 

3zio  - z\ o - -Zq o - 3z0o  + 

l2s  _ 
2 <'10 

lzs 

9 "00 

1 r 

2 2,31 

4 

= 

2Zoi  ~ 22o°  + ^zoi  _ 2oi  “ 

3 s 

2 "00 

3 zoo 

-5*10 

4i 

= 

— 3zoi  — ^oi  — z\\  + 3zn  — 

1 s 

22n  + 

1 5 

2 "01 

+ z00 

— 3zio  + z[0  + 3zoo  — 

, 1 s 
+ 22oo 

A 

c20 

= 

2Z\\  + 22i°  + 3201  + 2°J  _ 

3zn  + 

r'n 

izs 
- 2 2oi 

A 

c02 

= 

1 r 1 r o 3 5 

2Zn  ~ 2Zl°  ~ 3211  o2n 

i i zs 
+ 9 "01 

+ 3zio  + Zj0 
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c?2  — ^z\\  + 22l°  2ij  2Z°l  92°°  ^2qi 

+ ~oi  + 2oo  + 3^00  — 2 lo  — 3^io  + 3zn 

c30  = ~22”  ~ 92*0  ~ z°l  ~ 22°]  + 211  ~ 92o°  - 200  + 210 

C21  — — 3201  — 201  — Z\\  + 3211  — -2]]  + -2oi 

+ 2oo  _ 32io  + 2io  + 32oo  - 2 2io  + 22o° 

c03  = _2oi  - 211  + 2211  + 22®1  + 210  + 200  + 92i°  + 9 2oo 

c30  = 200  + 22oo  ~ 2n  + 2^11  “ 210  + 22i°  + 201  + 22®1 

c12  = ~ 9211  + 22*0  ~ 2^  ~~  ‘>2°1  + 920°  ^2°' 

+ 2oi  + 2oo  + 32oo  ~ 2io  — 3zio  + 3^11 

c03  = — 200  — ~2oo  + 2n  — 92ii  — 2io 

1 , , 1 s 

- 22l°  + 201  ~ 9 201 

6 SIbson  Element  Shape  Functions  on  a Unit  Cell 

In  this  section  we  introduce  the  coefficients  given  in  (95)  into  the  appropriate  polynomials  (1)  or  (2)  given 
in  Section  2.  We  will  then  rewrite  the  polynomials  in  the  form  of  shape  functions  where  the  coefficients  of 
the  twelve  corner  data  items  are  polynomials,  whose  order  does  not  exceed  three.  This  is  done  by  simply 
rearranging  the  spatial  and  corner  data  in  the  polynomials  (1)  and  (2)  after  substituting  (95).  The  Sibson 
element  shape  functions  for  each  of  the  four  triangles  are  given  by 

Zo(r,  s)  = (1  - 3s2  - 3r2  + s3  + 2 r3  + 3 rs2)  z00  + ^r3  + ^rs2  - 2 r2  - ^s2  + r'j  Zq0 

+ ^rs2  + ^s3  - rs  - ^s2  + s^j  2q0  + (-2r3  - 3 rs2  + s3  + 3r2)  z10 

+ (r3  + \rs2  - r2^  zr10  + ^-rs2  + ^s3  + rs  - ^s2^  2*0  (96) 

+ (— 3rs2  - s3  + 3s2)  20i  + ^-|rs2  + 2oi  + Q«3  - s2  + rs2^  zsQl 
+ (3 rs2  - s3)  zu  + 2n  + (~rs2  + ^3)  2h 

2l (r,s)  = (-(1  - r)3  + 3(1  - r)2(l  - s))  200  + f-|(l  - rf  + (1  - r)2(l  - s)^j  Zq0 

+ Q(1  - r)2(l  - s))  2q0  + ((1  - rf  - 3(1  - r)2(l  - s)  - 2(1  - s)3  + 3(1  - s)2)  210 

+ f”(l  ~ r)3  + (!  “ r)2(!  ~s)  + ^(1  - r ? - (1  ~ *0(1  - «))  2io 
+ f “(i  - r)2(l  -«)-(!-  s)3  + (1  - «)2^  2lo 

+ (-(1  - r)3  - 3(1  - r)2(l  - s)  + 3(1  - r)2)  20i  (97) 
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+ - r)2  - (1  - r)2(  1 - s)  + (1  - r)2^  zT0l  + Q(1  - r)2(l  - s)  - ^(1  - r)2^  zsQl 

+ ((1  - rf  + 3(1  - r)2(l  - s)  + 2(1  - s)3  - 3(1  - r)2  - 3(1  - s)2  + 1)  zu 

+ ("U  - rf  - (1  - r)2(l  - s)  + |(1  - r)2  + (1  - r)(l  - s)  - (1  - r))  zru 

+ (~U  - r)2(l  - 5)  - (1  - s)3  + i(l  - r)2  + 2(1  - s)2  - (1  - s))  zsu 


z2(r,s)  = 


23  (r,s)  = 


(3(1  - r)(  1 - s)2  - (1  - s)3)  zoo  + Q(1  - r)(  1 - s)2)  zr00 

+ Q(!  - 0(1  - s)2  - i(l  - s)3)  zs00  + (-3(1  - r)(l  - s)2  - (1  - s)3  + 3(1  - s)2)  zl0 

+ Q(!  - r)(!  - s)2  - |(1  - s)2)  do  + ("(I  - r)(l  - s)2  - |(1  - s)3  + (1  - s)2j  4 

+ (-2(1  - r)3  - 3(1  - r)(l  - s)2  + (1  - s)3  + 3(1  - r)2)  z0i 

+ ^-(1  - r)3  - |(1  - r)(l  - s)2  + (1  - r)2^j  zT01 

+ ^(1  - r)(l  - S)2  - ^(1  - S )2  - (1  - r)(l  - S)  + ^(1  - sf^j  2q! 

+ (2(1  - r)3  + 3(1  - r)(l  - s)2  + (1  - s)3  - 3(1  - r)2  - 3(1  - s)2  + l)  zu 
+ ^-(1  - r)3  - ^(1  - »0(1  - s)2  + 2(1  - r)2  + i(l  - s )2  - (1  - r)\  zru 

+ ^-(1  - r)(l  - s)2  ~ ^(1  ~ sf  + (1  - r)(l  - s)  + |(1  - s)2  - (1  - s)^j  zfj 

(r3  + 3r2s  + 2s3  - 3r2  - 3s2  + l)  z0o  + ^r3  + r2,s  ~ f 7-2  ~ rs  + 2oo 
+ ^r2s  + 7-3  - - 2s2  + 2qo  + (~r3  _ 3r2s  + 3r2)  z10 


-f  ( ^r3  + r2s  ~r2)z\ 0 + ( -^r2s  + V 


‘10 


+ (r3  - 3r2s  - 2s3  + 3s2)  z0 1 + |^r3  _ r2s  _ ^r2  + 2oi 
+ Qr2s  + fi3  - *2)  2oi  + (~r3  + 3r?s)  zn  + Qr3  - r2s^j  zru  + ^r2s^ 


-li 


(98) 


(99) 


7 Shape  Functions  on  a Scaled  Grid 

In  the  previous  sections  we  constructed  the  Sibson  element  shape  functions  on  a unit  square.  In  this  section 
we  will  extend  this  construction  to  rectangular  grids  with  uniform  spacings  along  the  x and  y axes.  The 
Sibson  split  for  this  grid  spacing  is  given  in  Figure  2.  The  spacings,  although  uniform  along  each  axis,  need 
not  be  the  same  for  the  x and  y axes.  In  particular,  let  (x,,Xi+ 1)  x (y^,  yJ+1)  be  a subrectangle  in  the  overall 
mesh.  Since  mesh  spacings  need  not  be  uniform,  let  Ax,  = x,+J  — x*  and  A y0  = yJ+i  — y3  Now  we  define 

x - x, 

r ~ 

Ax, 
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y 


where  < x < xl+\  and 


and 


By  a similar  argument 


(x,-yJ+i) 


(x,+i-yJti) 


Figure  2:  Rectangular  Sibson  Element 


Vj  < Vj+ 1- 


Then 


„ = y-yj 
Ay3 


dx 

dr 

dy 

ds 


A Xi 
AV3 


zoo  — 


-oo 


dz 

dr 

dz 

ds 


dz  dx 

dz  dy 

^(Xi,y,)Ts=  AyjZ,r 


7T 

z10 

= AxiZi+l,j 

z\0 

~ AVjZi+l,j 

•vr 

2 01 

= A XizfJ+1 

7S 

Z01 

= AVjzlj+ 1 

Z\1 

= A Xizf+lj+l 

211 

= A yjZ?+1j+1 
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(100) 


(101) 


(102) 


(103) 


Also 


-00 

— 

zl,J 

-1(1 

= 

-01 

= 

-1 1 

= 

-l  + lj-f 

With  this  notation  the  cubic  shape  functions  for  the  unit  element  can  be  written  in  terms  of 
(xl,Xi+i)  x (yj.yj  + 1),  where  r and  s are  defined  by  (100). 


)(x,  y)  — (l  — .is"  — 3r“  + s3  + 2r'!  + 3rs2)  + A.r,  (?,!  + — rs~ 


Zr 


o""  + r 


+ A ijj  ( rs2  4 -s3  - rs  - -s~  4-  s | r3;  4 ( — 2r3  — 3 rs~  4-  s3  + ,ir~ ) c,+  i,j 


4-  A.r,  I 7’ 3 4-  -rs2  - r" 


-i  + i,  j 


1 s . 1 , 


A tjj  I — r.s‘  4 —s'  4 rs  — 


V.v 

-i+i, j 


4-  ( — 3rs“  — s'  4-  3s~)  c,,j  + i 4 A.?',  ( — — rs“  + — s 


1 .2  , 1 A 


-l.J+1 


4- A//,  (-s3  - S"  4 rs”  j ~f,j  + i 4 (3rs‘  - s3)  + +i..j  + i 

+ A.r,  ( -rs- j ;,'r+1,J  + i 4 A.v;  (-rs~  + -s3j  A+i.j+i 

ZiU’.f/)  = (-(1  - r)3  + 3(1  - r)2(l  - s))  zl} 

+ A.r,  (-±(1  - v )3  4(1-  r)2(l  - s)\  ^ + A//,  (1(1  - r)2(l  - s)J  A, 

+ ((1  - r)3  - 3(1  - r)2(  1 - s)  - 2(1  - s)3  4-3(1  - s)2)  2I+lj 
4 A.r,  (-1(1  - r)3  4(1-  r)2(l  - s)  + 1(1  - r)2  - (1  - r)(l  - s)j  <+lj 

+ A y,  ( -1(1  - r)2(  1 - s)  - (1  - s)3  4(1-  s)2j  <+lj 
4-  (-(1  - r)3  — 3(1  -r)2(l  -s)  + 3(l  - rf)  zlJ  + i 
+ A.r,  (-1(1  - r)2  - (1  - rf(  1 - s)  + ( 1 - r)2j  <J+I 

+ A//,  (l(l-r)2(l-s)-l(l-  r)2)<J+1 

+ ((1  - r)3  4-  3(1  - r)2(  1 - s)  + 2(1  - s)3  - 3(1  - rf  - 3(1  - s)2  4 1)  + +i,./+i 

-4  A.r,  (-1(1  - r)3  - (1  - r)2(l  - s)  + 1(1  r)2  + (1  - r)(l  - s)  (1  r)J  ;,!+ 

4 A/Zj  (-1(1  - r)2(l  - s)  - (1  - s)3  4 1(1  rf  4 2(1  - s)2  (1  - s)J  A+i,.,  + i 

-2^,?/)  - (3d  - r)(  1 - s)2  - (1  - s)3)  + A.r,  (1(1  - r)(l  s)2j  zJ,) 

+ A.Vj  (1(1  - r)(l  -s)2  - 1(1  -s)3) 
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(104) 

(•r.//)  f- 

(105) 


(100) 

i j+i 


23  (x,y)  = 


+ (— 3(1  — r)(l  — s)2  — (1  — s)3  + 3(1  - s)2)  zi+ij 

+ Q(2  “ r)(X  ~ s)2  ~\(l~  5)2)  zi+hj 

+ AV]  (-i1  ~ r)U  ~ s)2  ~ - sf  + (!  “ 5)2)  zi+i,j 

+ (-2(1  - r)3  - 3(1  - r)(l  - s)2  + (1  - s )3  + 3(1  - r)2)  zitj+1 

+ Ax,  ^-(1  - r)3  - |(1  - r)(  1 - s)2  + (1  - r)2^  zfJ+1 

+ AVj  ^(1  ~ »■)(!  - s)2  - ^(1  - s)2  - (1  - r)(l  - s)  + ^(1  - s)2^  z?j+ j 
+ (2(1  — r)3  + 3(l  — r)(l  — s )2  + (1  — s)3  — 3(1  — r)2  — 3(1  — s)2  + l)  z,+ i j+i 
+ Axi  ^-(1  - r)3  - |(1  - r)(l  - s)2  + 2(1  - r)2  + ^(1  - s)2  - (1  - r)^  z?+1J+1 

+ A2/?  (-(!  - r)(l  - s)2  ~ ^(!  - s)3  + (1  - r)(l  - s)  + ^(1  - s)2  - (1  - s)^j  z\_ 

(r3  + 3 r2s  + 2s3  — 3r2  — 3s2  + l)  z,j  + Ax,  ^r3  + 

+A yj  ( \r2s  + r3  - ~r2  - 2 s2  + s^  zfj  + (-r3  - 3 r2s  + 3r2)  zi+ ij 


(107) 


j+ij+i 


r2s  _ |r2  _ rs  + r j ^ 


+Ax,  ( ^r3  + r2s  — r2 


) Zi+hj  + A%  + ^2)  Zi+l,J 

+ (r3  — 3r2s  — 2s3  + 3s2)  zlJ+ 1 + Ax,  ^-r3  — r2s  ~ 2 + rs^  2f, 

+Ayj  Qr2s  + s3  - s2^  zyi  j+l  + (-r3  + 3r2s)  zi+1J+1 


(108) 


+Ax,  ( ^r3  - r2s  ) z3 


i+ij+i  + ( ~2r2s  ) zi+i,j+i 


8 Derivatives 

Since  Lavery  [4]  uses  the  first  and  second  derivatives  of  the  Sibson  element  splines  in  his  interpolation  and 
approximation  procedures,  we  will  include  both  sets  of  derivatives  in  this  section  for  completeness.  Again 
we  will  use  the  notation 


and  Ax  = %i+\  — x,,  Ay  = yJ+ 1 - yj.  Then 


X — Xi 

r — 

A Xj 

y - Vo 

AVj 

dr 

l 

dx 

Ax, 

ds 

1 

dy 

AVj 

(109) 


(110) 
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8.1  First  Derivatives 

From  the  cubic  splines  for  each  of  the  triangles  in  a cell,  given  by  (105),  (106),  (107),  and  (108),  we  can  apply 
the  chain  rule  and  (110)  to  compute  the  first  partial  derivatives.  Direct  computation  gives  the  following. 


dz0 

dx 


(x,y) 


dz°i  \ 

7^{x'y) 


dzi 

ox 


dz\ 

dy 


(x,y) 


„2\  yy 


_L  (_6r  + 6r2  + 3s2)  Zij  + (l  - 4r  + 3r2  + ^ zfj  + ^ (s  + s'2) 

+ ^6r  ~~  6r2  ~ 3s2^  Zi+lj  + (_2r  + 3r2  + 2 s2)  Z'+lJ  + Ax~  ~~  S ^ Z*+1'J 

+ (~3s2)  ZlJ+' + (4s2)  2b+‘ + ty (*2) 


il  („2\  y y 

~i,j  + 1 


(111) 


+ Ax~  (3S^  Z<+1J+1  + ( ~os 


S+l,j  + l 


2\ 


'i+l.j+l 


/\  y - / 3 \ 

— — (-6s  + 6rs  + 3s2)  + — — {-s  + rs)  zfj  + ( 1 - r - 3s  + 2rs  + -s2  ) zXJ 

+ ^ (-6rs  + 3fi2)  2i+i.j  + 7^7  (rs)  2?+u  + (r  ~ s - 2rs  + ^2)  2?+ij 

+ -^—  (6s  - 6rs  - 3s“)  Zj^+i  + (s  - rs)  zfJ+)  + 2s  + 27's  + -s2  J zfJ+1  (112) 

+ ^6r5  ~~  ^ 2l+1’J  + 1 + (_rs)  ^+bj+i  + (~2rs  + ^2)  2f+i,J+i 


(3(1  - r)2  - 6(1  - r)(l  - s))  ztJ  + ^3(1  - r)2  - 2(1  - r)(l  - s)^j  zxXJ 

+ _ 7’^1  _ + ^7  (_3(1  ~ r)2  + 6(1  - r)(!  - s))  zi+  U 

+ f~(  1 - r)  + |(1  - r)2  + (1  - s)  - 2(1  - r)(l  - s)^  zf+lj  + ((1  - r)(l  - s))  zy+lj 

+ -^T  (-6(1  - r)  +3(1  - r)2  +6(1  - r)(l  - s))  zliJ  + 1 (113) 

+ f — 2(1  - r)  + 3(1  - r)2  + 2(1  - r)(  1 - s)^  zfJ+1  + ^ ((1  - r)  - (1  - r)(l  - s))  zyJ+l 

+ ~Kx  _ _ 3^  _ r ^ “ r)0  “ s))  2‘+bj  + i 

+ f 1 - 3(1  - r)  + |(1  - r)2  - (1  - s)  + 2(1  - r)(l  - s)^j  z*+hj+1 

+ |j  (-(1  - ^ + (1  - r)(l  - s))  zyl+l  3+l 


_L(_3(i_r)2)2,J  + ^(_(i-+KJ  + (4(i-+) 

+ -T—  (3(1  - r)2  - 6(1  - s)  + 6(1  - s)2)  2i+1  J + 


/\  ’77 

— (3(1  - r)2  - 6(1  - s)  + 6(1  - s)2)  2i+ij  + — 1 ((1  - r)  - (1  - r)2) 


’i+  l,j 


+ ( ^(1  — r)2  — 2(1  — s)  + 3(1  — s)‘ 


^i^  + ^;(3{1~r)2)^+1 
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(114) 


+f| ((1  “ r)2)  *+' + (“5(1  “ r)2)  <*+■ 

1 /\  7*  • 

+ Ay~  (_3^  ~ r)2  + 6(X  _ s)  “ 6(J  _ s)2)  zi+ U+i  + ^ (~(1  “ r)  + (1  - r)2)  *?+ij+i 

+ (i  + ^(!  - r)2  - 4(!  - ^)  + 3(1  - s)2)  Zi+ij+i 


dz2 

dx 


(x,y) 


("3<1  ~ r)2)  *«  + H°  ’ r)2)  Z'>  + <-(1  ' r)2>  2« 

+ (3(1  “ r)2)  2l+1J  + (~5(1  “ r)2)  ^ + If  <(1  ^ r)2>  4+w 


+ ^T  (_6(1  - r)  + 6(1  - r)2  + 3(1  - s)2)  zlJ+ 1 + (-2(1  - r)  + 3(1  - r)2  + -(1  - s)2^  z?j+i 
+^T  ((!  - r)  - (1  - r)2)  zyJ+l  + ■£-  (6(1  - r)  - 6(1  - r)2  - 3(1  - s)2)  zl+ 1J+1  (115) 

+ (l  - 4(1  - r)  + 3(1  - r)2  + \{\  - s)2)  *?+1J+1  + (-(1  - r)  + (1  - r)2)  zf+1>J+1 


^(®.y)  = (-6(l-r)(l-s)  + 3(l-r)2)^  + ^(-(l-r)(l-s))2* 

+ ^ 2(1  - r)(l  - s)  + |(1  - s)2^  zyl3  + (-6(1  - s)  + 6(1  - r)(l  - s)  + 3(1  - s )2)  -zi+i,j 

((1  - a)  - (1  - r)(l  - s))  z*+1J  + (-2(1  - a)  + 2(1  - r)(l  - 8)  + |(1  - s)2)  zy+lj 

1 /\  y ■ 

+ •—  (6(1  - r)(l  - 5)  - 3(1  - s)2)  Zij+i  + — 1 ((1  - r)(l  - s))  zfJ+1  (116) 

lAyj  Liyj 

+ ((1  - r)  - (1  - s)  - 2(1  - r)(l  - s)  + |(1  - s)2^j  zf>j+1 

+ ^-  (+6(1  -r)-  6(1  - r)(l  - s)  - 3(1  - s)2)  zi+ ltj+1 
At 

+ A^(_(1~s)  + (1"r)(1's))^+1 
+ (l  - (1  - r)  - 3(1  - s)  + (1  - r)(l  -s)  + |(1  - s)2^  zf+1J+1 


(-6r  + 3r2  + 67-5)  zl0  + (l  — 3r  + ^ r 2 — s + 2 rs^j  2+  + (_r  4.  rs)  zj 

+~Kx~  (6r  _ 3r2  ~ 6rs) Zi+hj  + (-2r  + \r 2 + 2rs)  Z*+1J 


+K7 (r  “ rs)  z'+l'i  + 7TT  (3r2  ~ 6rs)  ^+1 


i J + 1 


+ (-r  + ^r2  + s - 2rs^j  zfJ+1  + ^ M < 

+ ^T  ^_3r2  + 6rS^  2j+1’J+1  + (|r2  ~ 2rS)  *<+U+i  + (_rS)  a:<+iJ+1 


(117) 
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dz3  \ 
dy 


J-(^-6s  + 6s>)zv  + ^{-r  + r’)z’, 

+ (l  + \rl  - is  + 3s2)  z\s  + (-3r2)  *+,„■  + ^ (r2)  *f+1J 

+ (_5r2)  2*,+1J  + + 6s  " Cs2)  2'J+' 

+tj- (r  “ r2)  2*?2+i  + (ir2  “ 2s  + 3s2)  <2+>  + (3r2)  2,+IJ+l 


(118) 


+ 


A^t  / ,2\  * 


Aj/j 


't+l.j  + 1 


+ ( 2* 


i+l.j  + l 


8.2  Second  Derivatives 

In  this  section  only  those  terms  with  nonzero  coefficients  will  be  written  out.  Thus  all  coefficients  of  corner 
data  items  that  are  zero  will  not  be  included  in  the  second  derivative  expansions. 


d2z0 
dx 2 


{x,y) 


d2  z() 
dydx 


{x,y) 


d2z0 

dy2 


(x,y) 


^(-6  + 120^  + — (-4  + 6 r)z* 

l L 


+ 


^2  (6  - 12r)  2i+l,J  + (~2  + 6r)  2i+lj 


Ax2 

1 


zzz^ (fo)  212  + ^ (s)  2« + zz <_1  + 2s)  2« 

+ (~6s)  2,+1J  + ^ (s)  2?+1'2  + zz (I  * 2s)  2 

+zzz^  (~6s)  2,j+i  + zjj (_s)  2fj+I 

+zz (2s)  <*' + a^~  (6s)  2i+i-2+i 

+^(“s)2f+1'2+1  + AZ(“2S,2?J+‘ 


~~~2  ( 6 T 6r  T 6s)  zl3  + 

^y)  : a 2/2 


Al‘(-l  + r)4 


Ay, 

+ — (-3  + 2r  + 3s)z?-  + ^2  (_6r  + 6s)^+ib 

/\  /T*  1 

+^M2^  + A^;("1“2r+35)2^ 

1 ,/\  T • 

+ A^2  (6  _ 6r  “ 6s)  2^+1  + (1  “ r)  ^J+1 


X—  (-2  + 2r  + 3s)  zf  +1  + ~r~2  (6r  ~ 6s)  zi+hj+1 
A y3  Ay; 

At  I 

+ A^f  {"r)  + A^  (“2r  + 3S)  Wi 


(119) 


(120) 


(121) 
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d22i , v 


^2  (“6(1  ~r)+  6(1  - s))  Zij  + — ^ (-3(1  - r)  + 2(1  - 5))  z. 

+^j| ((1  _ s))  zv  + zd2 (6(1 " r)  ~ 6(1  ~ s))  Zi+1J 


ij 


“ (2  — 3(1  — r}  — 2(1  — s))  2*j+i 


+ ~^T  (!  “ 3(1  - r)  + 2(1  - s))  zf+]  j + ^ (-(1  - s))  2i+i, 

+ a^2  (6  ~ 6(1  - r)  - 6(1  - s))  2*J+1  + ^7 
+ Aj^  (_1  + ^ ~ ^j+i  + ^2  (“6  + 6(1  - r)  + 6(1  - s))  zi+ ij+i 
+ ^7  (3  — 3(1  — r)  — 2(1  — s))  zi+1  J+1  + (1  ~ (1  — s))  2f+i,j+i 


(122) 


^21 

dydx 


[x,y)  = 


’ (6(l-S))^  + ^-(2(l-S))4  + ^-((l-r))4 


+ 4;(~1+2<1~r))2f+,j  + a7 


Axt&yj 

+a757(-6<1‘s),*'»‘j  + a7 

+^J77  (-of1  - r»  2u+>  + ^ f-^1  - r»  2?++>  + 57  <(1  - r))  2«+<  (123) 

i-  (1  - 2(1  - r))  zf+JJ+I  + (-(!  - ’■>Ki+. 


AxiAyj 

1 

AxiAyj 


+^777<6(1“r>)2'+,J+‘  + 4 


a2*, 

dy2 


d2z2 
dx 2 


(*.y)  = 


Ay: 


-^j(6-12(l-s))2i+1|j^-(2-6(l-s))2|'+lfi 

“ (4  — 6(1  — s))  zf+1  j+1 


—2  (-6  + 12(1  - s))  2i+l,j  + l — — 


(124) 


(^,y)  = 


— (Il9) 
3u&rv 


^2  (6  - 12(1  - r))  2-j+l^“  <2  - 6(1  - r))  <J  + , 

^2  (-6+12(1  - r))j,+iJ+1J-  (4  -6(1  -r))zf+li+1  (125) 

7777  (6<I  “ s))  2,2  + 77 ((1  ‘ s))  25  + a7  (2(1  “ s))  2i 

^-(-2(1  - a))  zf+hJ  (126) 


— yj  — —i 

+A^  (~6(1 " S,)2,+1J  + 77  ((1  “ s,)  2f+1'2  + 57 

l 1 


922:2 
dy 2 


(*.3/) 


+7777 (~6<1  - s))  2,J+' + 77 (_(1  “ s))  *«+>  + 57 (_1  + 2(1  “ s,)  *’j+i 
+7777 (6(1  “ s))2,+lj+1  + 57  (-(1  _ s))  2'+u+‘ + 77 (1  “ 2(1 ~ s))2^+> 

1 At 

^77  (6(1  - r)  - 6(1  - s))  2ij  + -^4  ((1  - r))  2? 

— (2(1  - r)  - 3(1  - s))zf  + ~^r~2  (6  - 6(1  - r)  - 6(1  - s))  zi+1 
^Vj  J Ay * 


+ 


Ayj 

A 7*  1 

+ ~Aip  ^_1  + ^ ~ Z'+1’J  + A y~  ^ ~ ~ ^ _ ~ Z'+hJ 


,3 
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(127) 


<922,3 
dx 2 


<y) 


d2z3 

dydx 


(x,y) 


d2zs 

dy 2 


(*>y) 


+ 


1 /\  7‘ 

a^2  (-6(J  - r)  + 6(1  - s))  ^ (-(!  - 0)  <7+1 


^2  (-6  + 6(1  - r)  + 6(1  - s))  zi+  ij+i 


+ A^  (1  + 2(1  “ T)  ~ 3(1  “ S))  2^+J  + AyJ 

Ar  1 

1 (1  - (1  - r))z?+1J+1  + — (3-2(1  - r)  - 3(1  - s))z?+lj+1 


+ 


Ay2 


■ s)  z. 


ij 


^2  (-6  + 6r  + 6s)20  + — (-3  + 3r  + 2s)<  + (-1  + 

+ ^2  (6  - 6r  - 6s)  zi+hj  + (-2  + 3r  + 2s)  zf+lJ  + ^ ((1  - s))  zf+1 

+ ^2  (6r  - 6s)  2 6*,j+i  +Z^(-1+3r~  2<  <?+i  + « <+i 

+ ^2  (“6r  + 6<  2i+U+l  +X^(3r~  2<  <+l,J  + l + (-<  <+l  J + l 


_L_(6r)2lJ  + zL(-i  + 2,)z^  + zL(,)^ 


Ax,  Ay., 

L^z^’  + s: 


+ 


+ 


+ A^  <-6r,2'+1'J  + (2r)  + Axi  + 

1 (~6,')  2-'+'  + a^ (1  - 2r)  2«+'  + AS"  (r)  2«+' 

(6r)  2,+1J+i  + ^ (-2r)  zf+I  J+1  + — (->•)  *fj+i 
^ (-6  + 12s)  *<j  + A-  (-4  + 6s)  z” 

+ Ay1  ^ ~~  12<  2j,j+1  + ~Ay~]  ^ ~ 2 + 6<  <d+1 


AxlAyJ 

1 

Ax,  A yj 


(128) 


(129) 


(130) 
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